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The electronic structure of a LaNiC>3 bilayer grown along the [111] direction and confined between insulating 
layers of LaA103 is theoretically investigated using a combination of first principle calculations and effective 
multi-orbital lattice models. The LDA band structure is well reproduced by a tight-binding model for the Ni-e s 
orbitals defined on the buckled honeycomb lattice. We highlight peculiar properties of this model which in- 
clude almost flat bands as well as linear and quadratic band crossing points. The effect of local correlations is 
discussed within the LDA+[7 scheme and within the Hartree-Fock approximation for interacting multi-orbital 
lattice models. Over a wide range of interaction parameters we find that a ferromagnetic phase is energetically 
favored. We discuss the possibility of additional orbital order which could stabilize a spontaneous Chern insu- 
lator with chiral edge modes or a staggered orbital phase with a \/3 x \/~3 reconstruction of the unit cell. By 
studying an interacting nickel-oxygen lattice model we find that the stability of these orbitally ordered phases 
also depends on the value of the charge-transfer energy. Controlling the charge-transfer energy might therefore 
be an important step towards engineering exotic electronic phases in certain classes of oxide heterostructures. 



I. INTRODUCTION 

Design, growth and characterization of artificial structures 
of complex oxide materials has become an important field of 
research in modern condensed matter physics. 1 ' 2 The wide 
range of available materials with a multitude of physical 
properties carry great potential to engineer novel electronic 
devices 3 with desired functionalities, such as superconduc- 
tivity, ferromagnetism, ferroelectricity or multiferroic behav- 
ior. There are various physical mechanisms which mod- 
ify the electronic properties at interfaces and in superlattices 
as opposed to the bulk systems which can lead to atomic, 
electronic 4,5 and orbital reconstructions. 6 

A promising strategy to tune the electronic properties at 
interfaces and in superlattices of complex oxides is to grow 
the heterostructure along uncommon crystallographic direc- 
tions, such as along the [111] direction of a cubic perovskite. 7 
One interesting feature of such structures is the possibility to 
create a large variety of effectively two-dimensional lattices 
by sandwiching a precise number of layers of an "active ma- 
terial" between a good insulator. For example, a (111) bi- 
layer of a cubic lattice forms a (buckled) honeycomb lattice 
while three adjacent layers form a three-dimensional version 
of the dice lattice. 7-10 Realizing these lattice geometries is par- 
ticularly interesting in view of engineering electronic phases 
with topologically non-trivial band structures such as Chern 
insulators 11 or time-reversal invariant topological insulators 12 
with their associated protected surface modes. 7 ' ' 10 

In this article, we discuss the electronic structure of a (1 1 1) 
bilayer of LaNiC»3 sandwiched between sufficiently large re- 
gions of LaArC<3. LaAlC<3 is a wide band-gap insulator 13 
with an experimentally determined gap of E g = 5.6 eV while 
LaNi03 is a paramagnetic metal in bulk. The formal valence 
of nickel is Ni 3+ which corresponds to a partially filled 3d- 
shell with an electronic configuration t^e^, i.e. one electron 
in the doubly degenerate e g manifold. The insulating LaAlC»3 
layers insure that the electronic degrees near the Fermi en- 
ergy are effectively confined to the buckled honeycomb lattice 
formed by the Ni ions of the (111) bilayer. Figure 1 shows 
the actual supercell (LaNiC<3)2/(LaA103)io used in our den- 
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FIG. 1. The supercell of the superlattice (LaNi03)2/(LaA103)io con- 
sidered in this work. The growth direction is the [111] direction of 
the cubic perovskite. The primitive cell of the perovskite structure is 
also indicated. 



sity functional theory (DFT) calculations. 

Digital heterostructures between LaNiOs and LaAlC<3 have 
been considered previously but the focus was on systems 
grown along the [001] direction. 16-19 The interest in the [001]- 
grown heterostructures has been stimulated (at least in part) 
by the prospect to control the orbital character of the con- 
duction electrons by applying tensile or compressive strain. 
As opposed to bulk LaNi03 where no sign of orbital pref- 
erence is observed, the reduced symmetry in the [001] het- 
erostructure splits the energy between the d z i and the d x 2_ y 2 
orbitals. Theoretically, it has been suggested to apply tensile 
strain to selectively lower the energy of the d x 2_ y 2 -orbital. In 
this case, one could hope to imitate the relevant conditions 
for cuprate high-temperature superconductivity with nickelate 
analogs. 20-22 However, instead of a phase with uniformly oc- 
cupied d x 2_ y 2 orbitals, a charge-density wave is experimen- 
tally observed for tensile strain, 16- similar to the insulating 
bulk phases observed in other nickelates RNiC»3 with R=Pr, 
Nd, or Sm. 23 

In contrast to the [001] heterostructure, the e g -orbitals of 
the Ni ions remain doubly degenerate for the heterostructure 
grown along the [111] direction as long as the trigonal sym- 
metry is preserved. The non-interacting band structure of an 
e g -orbital model has a degeneracy point at k = where two 
bands touch quadratically. Such a situations can in princi- 
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pie give rise to interesting weak-coupling instabilities, includ- 
ing spin nematic and topological insulator phases. 10 ' 24 From 
the theoretical point of view, strong electronic correlations, 
charge-transfer physics and the complex interplay between 
the electronic and lattice degrees of freedom pose many chal- 
lenges in modeling the nickelates. Moreover, to the best of 
our knowledge, there are currently no experimental data avail- 
able for the (111) bilayer system. In this study, we use natu- 
ral idealizations of the real system which allow us to under- 
stand the electronic structure from simple multi-orbital lattice 
models. However, due to the lack of phenomenological in- 
put, the present work should be considered as a first step to- 
wards a more accurate theoretical description of these inter- 
esting structures. 

This paper is organized as follows. In Sec. II we present 
DFT results of the electronic structure for the supercell in 
Fig. 1. We show that the electronic degrees near the Fermi en- 
ergy are predominately given by the e g electrons with mixing 
of the nearby oxygen-p electrons. Effects of local correlations 
are discussed within the LSDA+C/ scheme 25 which predicts 
ferromagnetism. In Sec. Ill we show that the band structure is 
well reproduced by a tight-binding model for e g electrons on 
the buckled honeycomb lattice and highlight universal prop- 
erties arising from the particular geometry. The effect of local 
correlations is considered within this reduced model using a 
conventional mean-field approximation. In agreement with 
the LSDA+J7 calculation, we find that a ferromagnetic phase 
is stabilized over a wide range of interaction parameters. De- 
pending on the relative strength of the Hund's coupling J, we 
also identify spontaneous orbital order. In Sec. IV we intro- 
duce an effective model which includes both the Ni-e ff and 
oxygen-p states. This allows us to study how the charge- 
transfer physics 26 alters the stability of the orbitally ordered 
phases. 



H. DENSITY FUNCTIONAL THEORY 

We have studied the electronic structure of the 
(LaNi03)2/(LaA103)io supercell (see Fig. 1) using den- 
sity functional theory (DFT) 1415 within the local density 
approximation (LDA) 15 as well as the local spin density 
approximation (LSDA) 27 as implemented in the Vienna 
ab-initio simulation package (VASP). 28 We used the projector 
augmented wave pseudopotentials for all our calculations. 29 
A plane wave cut off energy of 600 eV and a 6 x 6 x 6 fc -point 
grid was chosen for integrating over the Brillouin zone. The 
lattice constant for the supercell was chosen as 3.82 A which 
corresponds to the experimental pseudocubic lattice constant 
of bulk LaA103. Atomic relaxation effects were not taken 
into account for our calculations. 

The LDA band gap for LaAlC>3 is E g = 3.3 eV (experi- 
mental value Eg = 5.6 eV). 13 As expected, this wide band 
gap leads to a strong confinement of the electronic degrees 
of freedom to the LaNiC>3 bilayer. The spectral weight near 
the Fermi energy is predominately comprised of states located 
at the Ni and their neighboring O ions. On the adjacent Al 
ion, we find very little spectral weight near the Fermi energy. 



These states on the Al are mainly associated with the s and p 
orbitals which weakly hybridize with the O-p states connect- 
ing the nearby Ni ion. We do not find a recognizable differ- 
ence in the total charge of an Al ion located next to the inter- 
face or far away, indicating that there is essentially no charge 
transfer from the Ni to the Al. By symmetry, the "inner O 
layer" which is sandwiched between the two Ni layers is dis- 
tinguished from the two "outer O layers" adjacent to only one 
of the two Ni layers. Figure 2 shows the orbital projected for a 
Ni atom as well as for an inner and outer O as obtained within 
the LDA. The LSDA calculations gave identical results. In 
Fig. 2 we set the zero of the energy to the Fermi level. The oc- 
tahedral crystal field splits the Ni 3d orbitals into the triply de- 
generate t2 g states and the doubly degenerate e g states which 
are well separated near the Fermi level. The e g DOS near the 
Fermi energy has a three peak structure with two prominent 
features near E « and E < 2 eV. As we show in Sec. Ill, 
these peaks appear due to weakly dispersing bands which are 
generically present in an e g model on the buckled honeycomb 
lattice. The orbital projected for the O atoms show that there 
is considerable mixing between the Ni-e g and the O-p states. 
Out of the three p orbitals it is the p a orbital pointing towards 
the Ni ion which hybridizes most strongly. Also note that the 
outer O mixes strongly with the lower and the inner O with 
the upper flat band. The LDA band structure along a high- 
symmetry path in the Brillouin zone is shown in Fig. 3(a). 
Near the Fermi level, there is a well-separated block consist- 
ing of four bands for each spin projection. Comparing with 
the results for the projected DOS in Fig. 2 we can associate 
these four bands with the two e g orbitals and the A-B sub- 
lattice structure of the (111) bilayer. We will discuss further 
characteristic properties of this band structure in more detail 
in Sec. III. 

In order to study correlation effects of the d electrons in 
Ni we also performed LSDA+C/ calculations where we em- 
ployed the simplified rotationally invariant approach as in- 
troduced by Dudarev et al.- 5 We choose an effective U e g 
value of 5.74 eV which was taken from Ref. 30 where the 
t/eff parameter has been computed in a self-consistent ap- 
proach following the scheme introduced in Ref. 31. It has 
been shown that this choice leads to accurate results for the 
ground state atomic structure. 30 The chosen U c g is also con- 
sistent with the reported experimental value in a recent work 
by Nohara et. al 32 who performed GW as well as LSDA+t/ 
calculations on LaNiOs. Upon fitting the calculated energy 
spectra of a cluster model to experimental X-ray photoelec- 
tron spectroscopy (XPS) and X-ray absorption spectroscopy 
(XAS) data 33 an effective U eS of 5.7 eV was obtained. 

Our LSDA+C/ calculation favors a ferromagnetic state with 
a magnetic moment of 1.12/ie on each Ni atom, where {1b 
is the Bohr magneton. This is consistent with the theoreti- 
cal value of 1/ie from LSDA+C/ reported in Ref. 30 for bulk 
LaNi0 3 . Figure 3(b) shows the spin-resolved band structure 
near the Fermi level. The e g bands for spin-| and spin-| com- 
pletely separate. The Fermi level is in the fully polarized 
majority band and is now shifted to the vicinity of the linear 
band-crossing point (Dirac point) located at K (and K'). Fig- 
ure 4 shows the projected and spin-resolved DOS. Note that 
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FIG. 2. (Color online.) LDA result for the orbitally projected of the 
Ni ions and its two neighboring O ions which are distinguished by 
symmetry. 



the ferromagnetic solution remains gapless and is fully po- 
larized such that only the majority electrons contribute to the 
DOS at the Fermi energy. In our LSDA+C/ calculation there 
is no sign of orbital ordering. 

Bulk LaNiOa is experimentally found to be paramag- 
netic although LSDA+C/ calculations predict a ferromagnetic 
phase. 30 For the bulk system, it has been argued that the poor 
treatment of the dynamical screening in the LSDA+C/ is re- 
sponsible for the disagreement with experiments. 30,32 From a 
phenomenological point of view, U e g = gives the best agree- 
ment with experimental data 30 and this choice has been used 
as starting point for GW calculations. 32 The ferromagnetic 
state found in our LSDA+C/ calculation (with U c h = 5.74eV) 
for the bilayer system should therefore be taken with care. 
On the one hand, one expects that the static treatment of 
the atomic interaction could be even more problematic in the 
quasi-two dimensional bilayer system. On the other hand, or- 
dered phases not present in the bulk have been experimentally 
observed in strained LaNi03 thin films. 19 Due to the lack of 
phenomenological input from experiments, we think that this 




FIG. 3. (a) The LDA band structure near the Fermi level along a high- 
symmetry path in the hexagonal Brillouin zone, (b) The LDA+C/ 
band structure along the same path as in (a). The solid (dashed) 
curves denote the energy bands of the majority (minority) spin pro- 
jection. The zero of energy is set to the Fermi level. 



uncertainty can not be resolved at the present moment. 



m. e a -ORBITAL MODEL FOR THE (111) BILAYER 

Let us now discuss the LDA band structure from the per- 
spective of a tight-binding model defined on the buckled hon- 
eycomb lattice. The tight-binding model helps to identify 
generic features of the band structure in this particular geom- 
etry. Based on the observations made in the previous section, 
we focus here on an orbital model for e g electrons alone. This 
is justified if the crystal field splitting is sufficiently large and 
the energy of the oxygen p-states is well separated from the 
energy of the Ni-e g states. As we will show in the follow- 
ing, this approach allows us to reproduce the LDA band struc- 
ture fairly well. Nevertheless, the e g model discussed in this 
section misses an important physical ingredient, namely the 
charge transfer between the Ni and O ions. In Sec. IV we 
therefore review the role of the charge transfer in the context 
of a more general model which also includes oxygen p states. 



A. Idealized tight-binding model 

In the simplest tight-binding model we consider only the 
dominant hopping between neighboring transition-metal ions 
which is mediated via an intermediate oxygen p-orbital. It is 
convenient to define the following orbitals 
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e u -d, d= (d x 2_ y 2,d 3z 2_ r 2), 



(1) 



where e x = (VS/2, -1/2), e y = (-- v/3/2, -1/2, ) and e z = 
(0,1) are the unit vectors in the directions of the nearest- 
neighbor bonds projected to the plane perpendicular to the 
[111] direction. The real orbitals d x 2 and d y 2 are rotations of 
the d z 2 orbital around the [111] axis by ±2ir/3 and are shown 
in Fig. 5(a). For oxygen mediated hopping in the direction e u , 
the Slater-Koster energy integrals 34 are only finite between 
neighboring d u 2 -orbitals. The tight-binding model takes the 
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FIG. 4. (Color online.) LDA+17 results for the orbitally projected 
in the ferromagnetic phase. Shown are the Ni 3d orbitals and the p 
orbitals for the two neighboring O ions which are distinguished by 
symmetry. Spin-| and spin-^ contributions are shown separately. 



simple form H t = Y,a- Hta where 

zeA u=x,y,z 



+ h 



■c). (2) 



The operators d 2 creates an electron with spin a in a 
d u 2 -orbital at site i and the summation runs over the sites of 
the A sublattice of the honeycomb lattice. Employing peri- 
odic boundary conditions, Eq. (2) is readily diagonalized in 
momentum space. 7 The band structure shows two dispers- 
ing bands which are reminiscent of the electronic structure 
of graphene, 



e 2 (k) = -t/2^3 + A k , e 3 (k) = t/2^3 + A k . 
Here, we have introduced 

A k = 2cos(V3k x ) + 4cos(-^-fcx) cos(— — ). 



(3) 



(4) 



The subscript of the wave-vectors refers to the two- 
dimensional (X, y)-coordinate system, see Fig. 5, and the 




FIG. 5. (a) Oxygen mediated hopping along the principle axis of 
the cube are between the d 3u 2_ r 2 orbitals with u = x,y,z. (b) A 
spatially localized eigenstate of the nearest-neighbor tight-binding 
model Eq. (2). The viewpoint is along [111], (c) Taking equal weight 
superpositions of the hexagonal states (i) allows one to construct lo- 
calized states with support on arbitrary contractible loops (ii). There 
are two more states which are linearly independent from the hexag- 
onal states and have support on a loop encircling the torus in one or 
the other direction (iii). 



unit of length is chosen as the length of the projection of the 
nearest-neighbor bond into the [111] plane. The two dispers- 
ing bands cross linearly at K + = f^giO) ar, d K- = ~K+ 
forming Dirac cones as in graphene. However, different from 
graphene, the dispersing bands are "sandwiched" by two flat 
bands with energy 



ei(fe) = -3t/2, e 4 (fc) = 3i/2. 



(5) 



The two peaks observed in the projected of the Ni e g orbitals 
in the LDA calculation, Fig. 2, can be related to these two fiat 
bands. The flat bands touch the dispersing bands in a single 
point at k = 0. Remarkably, the model Eq. (2) is formally 
equivalent to the planar p-band model on the honeycomb lat- 
tice introduced in Ref. 35. 36 However, the physical implemen- 
tation of Eq. (2) is very different in the present context: in- 
stead of p-orbitals, the elementary degrees of freedom are the 
e s -orbitals of the <i-electrons which hop on a buckled honey- 
comb lattice formed by a cubic (111) bilayer instead of the 
planar hexagonal lattice. 

We can use this formal analogy to the planar p-model to 
give a physical interpretation of the flat bands and the degen- 
eracy points. This is done by constructing spatially localized 
eigenstates of Eq. (2). ,37 Owing to the fact that hopping in 
the e u direction is only finite between d u 2 -orbitals, it is pos- 
sible to find spatially localized eigenfunctions which extend 
around a single hexagon located at R, see Fig. 5(b): 
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(6) 
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Here, v = ± labels a member of the bottom (-) and top (+) 
flat band manifold and 6j = (j - 1)tt/3. At each site j of 
the hexagon, the orbital wave function is orthogonal to the 
"outward-pointing" orbital eL u 2, which confines the electron 
to the hexagon. Assuming N unit cells, there are a total of N 
hexagon states for a given v. Taking equal weight superpo- 
sitions of different hexagonal states, one can construct eigen- 
states of Eq. (2) with support on arbitrary contractible loops, 
see Fig. 5(c). 

From the definition of the hexagonal states Eq. (6) it fol- 
lows that Y,r \^r) - 0. Therefore, only N — 1 states are lin- 
early independent. However, one can construct two additional 
eigenstates in real space which are linearly independent from 
the states in Eq. (6). They have support on non-contractible 
loops which encircle the torus in one or the other direction, 
as shown in Fig. 5(c). Therefore, the total number of linearly 
independent eigenstates with energy equal to the fiat band is 
N+l, It follows that there must be an isolated point where the 
flat band touches the dispersing band. 35,37 The corresponding 
k = Bloch states are given by 

\^O.x 2 -y 2 ) = —j=Ff ( Y \ d r,x 2 -y 2 ) ~ v Y \ d i,x 2 -y 2 ) J > 
V/V \ieA inB I 

VVV \ieA inB I 

with Af a normalization factor. This result from real-space 
topology is consistent with the band structure obtained by di- 
rect diagonalization. 

B. Generalized tight-binding model 

The idealized model Eq. (2) of the previous section only in- 
cludes the hopping t between the directed <i M 2 -orbitals. To get 
a satisfying fit to the LDA band structure, the tight-binding 
model has to be generalized by including additional hop- 
ping terms. Although these terms are small they neverthe- 
less modify the band structure notably. In particular, ad- 
ditional hopping amplitudes introduce a weak dispersion in 
the otherwise flat bands. It turns out that the two-fold de- 
generacy at k = is protected by the trigonal symmetry 
(D^d) of the bilayer system and therefore is robust against 
the inclusion of additional tight-binding parameters. In fact, 
the two pairs of degenerated eigenstates at k = span the 
vector spaces of the two-dimensional irreducible representa- 
tions of the D^d point group. This fact is interesting because 
if the Fermi surface coincides with the symmetry-protected 
quadratic band touching point, unusual weak coupling in- 
stabilities are possible. 10 ' 24 ' 38 ' 39 The symmetry -broken state 
generically is either a (lattice) nematic phase in which the 
rotational symmetry is broken or a topological phase with a 
gapped band structure characterized by a non-trivial topolog- 
ical invariant. 

In the following we refine the tight-binding model by in- 
cluding additional parameters consistent with the e g character 
of the orbitals and assuming that the trigonal symmetry of the 
lattice is preserved. 



1. Nearest-neighbor hopping 

The Slater-Koster parameters for hopping along the z- 
direction yield the matrix 



\0 t s ) 



(8) 



in the basis (d z 2,d x 2_ y 2). As before, t includes predomi- 
nantly the hopping via the intermediate oxygen while tg arises 
from the direct overlap and is small, see below. Assuming 
that the nearest-neighbor hopping in x and y directions are 
equivalent to the hopping along the z direction, we obtain the 
corresponding matrices by a rotation of the e g -orbitals around 
[111] by ±2n/3. The matrix for the rotation by 2n/3 is 



R 



As a result, we then find 



-1/2 n/3/2\ 
-V3/2 -1/2/ 



(9) 



(10) 



2. Second-neighbor hopping 

The Slater-Koster parameters for second-neighbor hopping 
via two intermediate oxygen atoms define the matrix 



t 



t'/2 
-V3A/2 



V3A/2 
-3t'/2 



(11) 



They take into account the lowest-order processes for the 
second-neighbor hopping. The off-diagonal entries propor- 
tional to A are allowed in the bilayer system discussed here 
(as opposed to a perfect cubic system) because the two pos- 
sible paths connecting second-neighbor transition-metal ions 
are not equivalent: they either involve "inner" or "outer" oxy- 
gens, see also Fig. 9. Note that t xy is not symmetric if A + 
which means that there is an associated direction for the hop- 
ping. We use the convention that t xy denotes the hopping of 
an electron along a second neighbor bond which is reached 
by first following the y-axis and then the x-axis of the cube. 
By rotating the orbitals, we also obtain the second-neighbor 
hopping along the other directions: 



tyz — R t X yR^ i zx — R ty Z R. 



(12) 



Including the above introduced hopping matrices, the gener- 
alized tight-binding model now takes the form 

H 0= YY Y (dl,Judcr,r+e u + h.C.) 
rzA £7 u=xyz 

+ YY Y {di,Ju,u+ld<T,r+e u -e u+1 + h.C.) (13) 
reA <r u=xyz 

+ YY Y (dl,riu,u+ld<r,r-e u +e u+1 +h.C.) . 
rzB <7 u=xyz 

Here, d a = (d z 2_ a ,d x 2_ y 2 a ) T is a vector in orbital space and 
the notation u + 1 refers to y if u = x with a cyclic extension 
to the other elements. 
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FIG. 6. (Color online). The solid line shows the LDA bands near 
the Fermi energy along a high-symmetry path T-K-M-r for an unre- 
laxed (LaNi03)2/(LaA103)io superlattice. Also shown are the tight- 
binding bands obtained from a least-square fit as detailed in the main 
text and Tab. I. Shown are the fits A (diamonds), B (squared) and C 
(circles) while fit D is essentially indistinguishable from C and not 
shown. The bands of A have been shifted by +0.2 eV for better visi- 
bility. 



fit 


UeV] 


*' [eV] 


A [eV] 


ts [eV] 


E F [eV] 


A 


0.603 











-0.701 


B 


0.600 


0.058 








-0.693 


C 


0.598 


0.062 


-0.023 





-0.693 


D 


0.598 


0.062 


-0.023 


-0.007 


-0.693 



TABLE I. Parameters obtained in different tight-binding fits to the e g 
LDA band structure of the 12 layer superlattice shown in Fig. 6. with 
an increasing number of adjustable parameters. 



C. Tight-binding fit of LDA band structure 

We used the refined tight-binding model Eq. (13) for a least- 
square fit of the LDA band structure of the e s -bands along the 
path T-K-M-r. The result is shown in Fig. 6 and the fitting 
parameters are listed in Tab. I. We have performed fits with 
an increasing number of adjustable hopping amplitudes t, t', 
A and tg. We find that the nearest-neighbor hopping is about 
t = 0.6 eV. The next biggest parameter is the second neigh- 
bor hopping for which we find t' « O.lt. This values are in 
agreement with previous findings. 8 ' 40 Even smaller is the off- 
diagonal term A in Eq. (11) which is specific to the bilayer 
system and we find A « -0.04t. However, the inclusion of 
this parameter markedly improves the fit, especially near the 
if -point, as can be seen in Fig. 6. Finally, we find that the in- 
clusion of the direct overlap tg is vanishingly small and does 
not lead to an essentially better fit. 

We have also performed fits whit a tight-binding model with 
completely general hopping matrices for t z [Eq. (8)] and t xy 
[Eq. (11)]. However, we found that the least-square fit is only 



marginally better indicating that the most dominant processes 
are included in the form of the matrices in Eqs. (8) and (11). 
Despite the good overall agreement, certain details of the LDA 
band structure (e.g. the low-lying part near M and along M-T) 
are not captured within the present model. This indicates that 
in order to describe all the details of the LDA band-structure 
in a e g tight-binding model, one has to consider further-range 
hopping parameters. An alternative way is to explicitly in- 
clude the oxygen-p states. This is done in Sec. IV where we 
discuss a Ni-O lattice model. 



D. Interaction effects in e g -model 

We now turn to an analysis of interaction effects within the 
effective description of the tight-binding model Eq. (13). Ow- 
ing to the localized character of the 3e?-orbitals, we include 
local interactions of the standard form 41 ' 42 

(u'-J) E 

r a a>(i,a 

+ U' E nratrirfU + J E 4 Q 

+ I E d l^ d rp\j ra l d rPl\ (14) 

We assume the following relations between the Slater- 
Kanamori interaction parameters: U' = U-2J and I = J. 
They are valid in free space and considered as approximately 
true in the solid state environment. The total multi-orbital 
Hubbard Hamiltonian for the e g electrons is given by 

H = H Q + H mt . (15) 

The interacting Hamiltonian Eq. (15) with only nearest- 
neighbor hopping has been studied previously within the 
Hartree-Fock approximation and the phase diagram has 
been worked out for various combinations of interaction 
parameters. 810 A particularly interesting result for intermedi- 
ate to strong interactions is the observation of a spontaneously 
generated quantum anomalous Hall (QAH) phase which is ac- 
companied by ordering in complex orbitals within a ferromag- 
netic (FM) phase. The resulting mean-field band structure is 
topologically non-trivial with a finite Chern number v = ±1 
displaying a spontaneous quantum Hall effect. 

Here, we generalize the previous Hartree-Fock studies in 
three ways: We study the effect of a finite second-neighbor 
hopping t', we allow for ordering with a tripled unit cell and 
we also consider larger values of the Hund coupling J. The 
resulting phase diagram for fixed interaction U + J = lOt as 
function of the ratio J/U and t' /t is shown in Fig. 7. Our new 
findings are summarized as follows: (i) For 0.1 < J/U < 0.2 
we find that a ferromagnetic phase (FM) with staggered orbital 
order (AFO) which triples the unit cell but preserves the three- 
fold rotation symmetry has lower energy than the previously 
found phase which breaks the rotation symmetry but pre- 
serves the original unit cell. 81 (ii) For large Hund coupling 
J/U > 0.4 a charge-density wave (CDW) is stabilized, (iii) 
For the considered interactions, there is very little dependence 
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FIG. 7. Hartree-Fock phase diagram obtained from the interacting 
e 9 -model for fixed interaction U + J = Wt as function of the ratio 
J/U and t'/t. The various phases are explained in the main text. 



1. FO/AFM 

For antiferromagnetically coupled layers, we found ferro- 
orbital (FO) order with 





W 



(18) 



Other equivalent orbital orders are obtained by rotating the 
real orbitals around [1 1 1] by ±2ir/3: 



(19) 



Here, the rotation matrix is given in Eq. (9). Hence, the 
FO/AFM order breaks the rotation symmetry while preserv- 
ing the translational symmetry. 



2. v^3 x %/3 AFO/FM 



on t' . As a side note, we remark that the phase-boundaries at 
weak interactions, as discussed in Ref. 10, are more strongly 
affected by the inclusion of the second-neighbor hopping t' . 
In particular, upon increasing t'/t, we observe a rapid shrink- 
ing (t'/t = 0.01 is enough) of the region where the sponta- 
neous topological insulator phase is energetically favored. 

In the following, we discuss the various phases shown in 
Fig. 7 in more detail. For the considered range of parame- 
ters, we find ferromagnetic order of the spin degrees within a 
Ni-layer. However, depending on the relative strength of the 
Hund exchange, the coupling between the two layers is anti- 
ferromagnetic (AFM) if J/U < 0.1 or ferromagnetic (FM) if 
J/U > 0.1. 

In addition to the magnetic order, we find various types of 
orbital order. The orbital order is conveniently discussed with 
the help of the orbital pseudospin-1/2 operator 



Ti = ~ E E dhaTapdipv. (16) 



The spontaneous development of long-range orbital order of 
the orbital pseudospin is indicated by a finite expectation 
value 



Mi = <^>- (17) 



Depending on the details of the interaction parameters, we 
find different types of orbital order. 



Increasing J/U we find a first-order transition to the ferro- 
magnetic phases. In this regime, spin-up and spin-down bands 
are separated by an energy (U + J)/2 and we find fully spin- 
polarized phases. Because of the large energy splitting be- 
tween t and I spins, the orbital order in this regime can be 
studied in a spinless model. Assuming polarization along the 
t-direction and neglecting the | -bands, the effective Hamilto- 
nian in the fully polarized ferromagnet is 

H FM = H 0tt +VY, n ia \n ibh (20) 

i 

with V = U'-J = U-3J. This model is formally equivalent to 
a "single-orbital" Hubbard model with "spin-dependent" hop- 
ping and Hubbard interaction V. The model Eq. (20) should 
be studied at half filling where the Fermi energy crosses the 
Dirac points for V = 0. 

The \/3 x \/3 AFO/FM phase corresponds to the intermedi- 
ate to large V limit of Eq. (20). It is characterized by a special 
type of staggered orbital order with a \/3 x \/3 reconstruction 
of the unit cell, see Fig. 8. In this phase, the orbital order pa- 
rameters at different sites are related by a unitary transforma- 
tion. For example, the order parameters at site 3 of Fig. 5(b) 
are related to site 1 in the following way: 

(*)-*(*) 

This corresponds to a rotation of the orbitals around the [111] 
axis by an angle of 2tt/3. A real space sketch of the orbital 
order in this phase is shown in Fig. 8. On the mean-field level, 
there is a sk-fold degeneracy of this phase. 

The findings of the mean-field theory can be compared with 
the strong coupling expansion of Eq. (20). At half filling, 
Eq. (20) reduces to an orbital pseudospin model of the form 

H K = KY, (rf rl ex + rfrf +ey + t*t^ ) (22) 
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FIG. 8. (Color online.) Sketch of the staggered orbital order phase 
with a tripled unit cell (\/3 x >/3 AFO/FM). The shown real orbitals 
indicate the majority occupation and the enlarged unit cell is shown 
as shaded area. 

where K = 2t 2 /V. In terms of the orbital pseudospin-1/2 
operator 2$, the operators entering Eq. (22) are given by 

rf = -l/2T t z -S/2T*, (23) 

rf = - 1 /2T? +V3/2Tf, (24) 

t? = T*. (25) 

Because K > 0, staggered orbital order is energetically fa- 
vored, consistent with the mean-field result. Nevertheless, the 
pseudospin interaction described by Eq. (22) is frustrated on 
the honeycomb lattice and to our knowledge the nature of the 
ground-state order hasn't been identified completely. Indeed, 
the study of the classical version of Eq. (20) in Ref. 43 shows 
that the ground-state is macroscopically degenerate. Further- 
more, it was shown that both thermal and quantum fluctua- 
tions lift this degeneracy. The \/3 x \/3 reconstruction identi- 
fied in our mean-field analysis is insofar a reasonable scenario 
because it maximizes the number of "resonating hexagons" 
favored by ring exchange processes. Such a behavior would 
be similar to the one observed in other frustrated models such 
as the hardcore dimer model on the hexagonal lattice 44 or the 
frustrated charge models on the kagome 45 ' 46 or checker board 
lattices 47 



3. QAH/FM 

In a narrow region of J/U, which corresponds to an inter- 
mediate effective interaction V in the spinless model Eq. (20), 
we find a phase which is characterized by the uniform orbital 
order along the y-direction, i.e. 

fly*0. (26) 

This means that the electrons predominantly occupy complex 
orbitals 

\d±id) = — {\d z *)± i\d x 2_ y 2 )). (27) 

Such orbital order opens a gap throughout the Brillouin zone. 
It has been pointed out previously 810 that the resulting insu- 
lator is a spontaneous quantum anomalous Hall (or Chern) in- 
sulator with topologically protected chiral edge modes. The 



mean-field band structure is characterized by a finite Chern 
number 48 n = ±1. 



4. Gapless FM 

Increasing the ratio J/U even further, we find a gapless 
ferromagnetic phase without long-range orbital order. This 
phase corresponds to the weak-coupling phase of the spinless 
model Eq. (20) and it reflects the well-known fact that the 
Dirac semi-metal is perturbatively stable against interactions. 
Therefore, orbital order is suppressed and fl = 0. Instead, the 
ferromagnetic solution with gapless single-particle excitations 
is stable. 



5. CDW/FM 

By further increasing the relative strength of the Hund cou- 
pling, we identify a second order phase transition to a charge- 
density wave (CDW) with charge disproportion between the 
top and bottom layer (i.e. between the A and B sublattices 
of the buckled honeycomb lattice). The tendency towards 
a CDW phase can be understood from the effective model 
Eq. (20) which shows that the effective interaction V becomes 
attractive if J > U/3. Our mean-field analysis reveals that an 
even larger critical ratio J/U ~ 0.4 is needed to stabilize the 
charge disproportion. 

A CDW with charge disproportion between the A and B 
sublattices has also been proposed for bulk nickelates as a 
mechanism to lift the orbital degeneracy which is different 
from the Jahn-Teller distortion. 49 Our findings that a relatively 
large Hund coupling is required is in agreement with a recent 
work based on a two-band model for bulk LaNi03. 40 



IV. NICKEL-OXYGEN MODEL 

In the previous section we discussed the band structure of 
the conduction electrons in terms of an effective multi-orbital 
Hubbard model for the e g states. In this model, the oxygens 
enter the description implicitly by mediating the hopping be- 
tween the Ni ions. The satisfactory fit of the LDA band struc- 
ture from the tight-binding model obtained in Sec. IIIC seems 
to justify this view -point. On the basis of this model, we have 
obtained the Hartree-Fock phase diagram which, depending 
on the interaction parameters, shows both orbitally ordered 
and disordered phases. While the ferromagnetic order pre- 
dicted by the Hartree-Fock for a large region of parameters is 
consistent with the LSDA+C/ result, we did not observe or- 
bital ordering in the DFT calculations. 

Below we reexamine the stability of the orbitally ordered 
phases in the context of a lattice model which explicitly in- 
cludes both the oxygen and nickel ions. Within the Hartree- 
Fock approximation, we study the dependence of the elec- 
tronic phases on the energy splitting between the Ni-d and O-p 
states and find that orbital order is more stable if the charge- 
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outer O 




FIG. 9. (Color online.) Sketch of the Ni-0 model for the (111) bi- 
layer. There are two layers of Ni ions denoted by Nil and M2 for 
which we keep the e g orbitals. There are an inner and two outer lay- 
ers of O ions. We keep the O-p orbitals which form a a-bond with 
their neighboring Ni-e 9 orbitals. The Al ions are treated as vacuum, 
so that there is no mixing with its neighboring oxygens. 
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transfer energy is increased. We discuss possible mechanisms 
leading to this result. 



A. Ni-O model for (111) bilayer 

In the following, we consider a model which treats both the 
Ni-d and O-p orbitals. 41 Instead of keeping the full Ni-cZ and 
O-p manifolds, we consider a simplified version with two e g - 
orbitals at the Ni sites and one p-orbital at the O sites, see 
Fig. 9. The Al ions are treated as vacuum. The p orbital is 
chosen such that it points along the directions of the cartesian 
coordinate system, i.e. the p-orbital which makes a a -bond 
with its neighboring d x 2_ y 2ld^, z 2_ r 2 orbitals. The resulting 
model has the form 

e piPi<rPi° + S e dd\ aa d ia(7 + H hyh + H p 

+H int + H DC . (28) 
The hybridization between the p and d electrons is given by 

#hyb= £ (Ol<W + h " C -)- (29) 

The hybridization V£- between the e g orbitals and the p or- 
bitals is parametrized by the Slater-Koster parameter (pda). 
H p describes the direct overlap between O-p orbitals for 
which we introduce the hopping parameter t pp . Finally, only 
the d-electrons are assumed to be correlated and H- m t is given 
by Eq. (14). Hdc accounts for the "double-counting" in the 
interacting model and will be discussed later. 

The unit cell of the bilayer system has two nickel and nine 
oxygen atoms arranged in layers along the [111] direction as 
Os/Ni/Oa/Ni/Oa. By symmetry, the inner oxygen layer, which 
is sandwiched between the Ni layers, has a different onsite- 
energy 4 than the outer oxygens e p °^ . 

B. LDA fit based on Ni-O model 

We find that there is an ambiguity in fitting the parame- 
ters of the Ni-O model to the LDA results. In particular, fit- 



FIG. 10. (Color online.) Comparison between the LDA and the Ni-O 
tight-binding (TB) model. In (a) we show the least-square fit of the 
e fl band structure along a high symmetry line. Panels (b)-(d) show 
the projected density of states obtained from the TB and the LDA 
calculations. The Ni-e s states are shown in (c) and the inner and 
outer O-p states in (b) and (d). Along with the total O-p DOS we 
also show the contribution from the O a orbital alone. 



ting different quantities accessible in the Ni-O model, such 
as the e g band structure, the orbitally projected DOS or any 
combination of these observables result in different optimized 
tight-binding parameters. Nevertheless, different fits which 
are in agreement with the overall structure of both the O and 
Ni states yield similar parameters. Our qualitative conclu- 
sions are therefore not affected by the details of the fitting 
procedure. Here, we show the result obtained by fixing the 
Slater-Koster parameters at (pda) = 1.8 eV and the hop- 
ping between oxygens at t pp = 0.7 eV. We choose to fit the 
e g band structure by optimizing the onsite energies for the 



and 



Ni states, e^, and for the inner and outer O states, e p 

e p . The fixed hopping parameters are similar to the ones 
used by other groups, 41 ' 50 and we find that this choice gives 
satisfactory agreement both for the e g bands as well as for 
the projected DOS. The optimized values are = -1.47 eV, 

4° = -4.74 eV and e ( p o) = -5.47 eV. The result of this fit is 
summarized in Fig. 10 where we show the e g bands and the 
projected DOS. 



C. Hartree-Fock approximation and double counting 

The Hartree-Fock decoupling of H lnt [Eq. (14)] yields a 
term 



(30) 



where (rid) is the averaged occupation of the Ni states and 
U = (3f7 - 5J)/4. In general, H c leads to a shift of the 
e?-orbital energy as compared to the O-p levels which would 
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FIG. 11. Dependence of the symmetry-broken phases on the shift in 
the d-electron level Aed as obtained from the interacting Ni-O model 
in the Hartree-Fock approximation. Shown are the average occupa- 
tion of the Ni-e 9 orbitals (rid), the uniform magnetization (m) of the 
ferromagnetic (FM) order as well as the expectation values of the or- 
bital pseudospin (T). Increasing Ae^ triggers a transition from the 
gapless FM phase to a gapped topological phase (QAH) and eventu- 
ally to a phase with staggered orbital order (AFO). Also shown is the 
bandwidth W (in eV) of the e g band obtained in the corresponding 
non-interacting Ni-0 model. 



modify the average occupations of the rf-orbitals when inter- 
actions are included. In the LSDA+C/ method, 51 one there- 
fore adds a double counting term which essentially compen- 
sate H c . It is argued that the DFT calculation already includes 
a local term similar to Eq. (30) and that the overall charge 
distribution (although not the orbital occupation) is obtained 
reliably from the DFT. The interacting Hamiltonian with the 
double-counting correction should therefore reproduce the 
DFT result in the absence of interaction-driven orbital order. 
This so-called double-counting problem appears also in the 
DFT+DMFT approach. 52 In the spirit of the LSDA+t/ we 
choose the double-counting term Hoc sucn that it compen- 
sates for the above discussed Hartree shift, Hoc = -He- The 
remaining interactions are treated in the conventional Hartree- 
Fock approximation and we solved the self-consistency equa- 
tions via iteration. 



Around Aed = 0, we find that a gapless ferromagnetic phase 
without orbital order is energetically favored. This is con- 
sistent with the LSDA+C/ results reported in Sec. II. Increas- 
ing Aed = leads to orbital ordering of complex orbitals and 
(T y ) + (QAH/FM). In this phase, the single-particle spec- 
trum is gapped but the band-structure is topologically non- 
trivial and characterized by a finite Chern number. By in- 
creasing Aed even further, we find that staggered orbital order 
(AFO) in real orbitals is favored and (T x ) + on each sublat- 
tice. 

There are two aspects of the charge-transfer physics which 
can help understanding these results. First, increasing the 
charge-transfer energy reduces the effective bandwidth of the 
e g bands. We illustrate this point in Fig. 1 1 by plotting the 
width W of the e g bands in the corresponding non-interacting 
Ni-O tight-binding model. The ratio U/W is increased by 
increasing Aed and correlation effects become more impor- 
tant. The observed orbital order can therefore be consid- 
ered as a result of increasing U /W in an effective model for 
the e g -electrons. 8 ' 10 A second aspect of the charge-transfer 
physics can be quantified by the average occupation of the Ni- 
e g states, (rid), which is also modified by the charge-transfer 
energy. For Ae^ = 0, the Ni-O model predicts (rid) w 2 which 
is very different from the ionic limit (rid) = 1 (correspond- 
ing to Ni 3+ ) and indicates that a large number of holes resides 
on the oxygen sites. In our LSDA+C/ calculation, the differ- 
ence from the ionic value seems even larger: examining the 
local density-matrix we find (n^DFT ~ 2.4. Previously, this 
effect was called "self-doping" of the oxygen p-band 53 ' 54 and 
is considered to capture aspects beyond modifying the ratio 
of U/W in an effective model for the e g electrons. Recently, 
it has been argued that (rid) is a fundamental quantity which 
crucially affects the interaction-driven metal-insulator transi- 
tion in the nickelates. 55 In Fig. 1 1, the reduction in (rid) g° es 
hand in hand with the reduction of W and therefore, these 
two aspects can not be separated in a clear way. Nevertheless, 
our results are consistent with the statement that the (gapped) 
orbitally ordered phases are more stable if fewer charges are 
transferred between the Ni and O. The charge transfer there- 
fore offers an explanation why the orbitally ordered phases are 
not observed in our LSDA+C/ calculation. 



V. CONCLUSIONS 



D. Charge transfer 

In order to see how the symmetry-broken phases are af- 
fected by the charge-transfer physics, 26 ' 53 ' 54 we performed 
Hartree-Fock calculations for different charge-transfer ener- 
gies using the double-counting term discussed in the previ- 
ous section. We fixed the tight-binding parameters at the 
values obtained from the LDA except for Ed which we var- 
ied according to Ed £d + Ae^. This effectively changes 
the energy splitting between the O and Ni states. Figure 1 1 
shows the dependence of the symmetry -broken phases on Ae^ 
for fixed interaction parameters U = 6 eV and J = 0.5 eV. 



In summary, we discussed the electronic structure of a 
LaNi03 (111) bilayer sandwiched between several layers of 
the insulator LaAlOs. Using a combination of first prin- 
ciple methods and effective multi-orbital lattice models we 
have studied both the non-interacting band structure as well 
as possible instabilities driven by local correlations among 
the Ni-d electrons. The non-interacting band structure is well 
reproduced by a generalized tight-binding model including 
also second-neighbor hopping amplitudes for e g states on the 
buckled honeycomb lattice. If local interactions among the 
Ni-ci electrons are included, both the Hartree-Fock and the 
LSDA+C/ predict ferromagnetic ordering over a wide range 
of parameters. We have extended previous Hartree-Fock 
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calculations 810 in various ways and, apart from the magnetic 
order, we have also identified possible orbital ordering which 
could give rise to interesting electronic phases including a 
spontaneous Chern insulator or a staggered orbitally ordered 
phase with a triplet unit cell. Furthermore, we identified a 
charge-density wave for large values of the Hund coupling. 

By considering an effective lattice model which includes 
both the O-p and Ni-d states, we studied the effect of mod- 
ifying the charge-transfer energy between the Ni and O 
states on the orbital ordering. Our results are consistent 
with the statement that orbital ordering is suppressed the 
more charge is transferred. This is in qualitative agree- 
ment with recent theoretical results on the orbital polariza- 
tion in [001] heterostructures 50 ' 56 and the phase boundaries 
for the paramagnetic metal-insulator transition in bulk nicke- 
lates and cuprates. 55 The charge-transfer energy was experi- 
mentally observed to change with n in (LaA103)3/(LaNi03)„ 



superlattices suggesting that it might be possible to control 
the degree of covalency in artificial structures. 19 Important is- 
sues such as the atomic relaxation in the bilayer geometry 
and its dependence on the electronic correlations as well as 
the treatment of the electronic correlations beyond the static 
mean-field approximations remain topics for future study. 
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